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Markov  Chain  Simulations  of  Binary  Matrices 

William  B.  Krebs 

Abstract:  We  consider  Markov  chains  to  simulate  graphs  with  a  fixed 
degree  sequence  and  binary  matrices  with  fixed  row  and  column  sums.  By 
means  of  a  combinatorial  construction,  we  bound  the  subdominant  eigen¬ 
values  of  the  chains.  Under  certain  additional  conditions,  we  show  that  the 
bounds  are  polynomial  functions  of  the  degree  sequences  and  the  row  and 
column  sums,  respectively. 


1.  Introduction: 

Let  M  be  a  given  m  x  n  matrix  whose  entries  are  Os  and  Is.  We  want  to  choose  a 
matrix  N  distributed  uniformly  over  those  m  x  n  0-1  matrices  with  the  same  row 
and  column  sums  as  M. 

Let  71,^2]  be  a  2  x  2  submatrix  of  M  taking  one  of  the  two  forms 


1 

o‘ 

0 

r 

0 

1 

or 

1 

0 

Exchange  Os  and  Is  in  M[ii^i2;ji^j2]i  suid  call  the  resulting  matrix  M'.  It  is  clear 
that  the  row  and  column  sums  M*  are  the  same  as  those  of  M.  the  operation  taking 
M  to  A/'  is  called  an  interchange. 

The  interchange  operation  defines  a  graph  structure  on  the  binary  matrices  whose 
row  and  column  sums  equal  those  of  M.  The  set  of  such  matrices  forms  our  vertex 
set;  if  N  and  are  vertices,  an  edge  joins  them  if  there  is  an  interchange  taking 
N  to  N'.  Call  this  graph  the  interchange  graph. 

One  way  01  si..''  lating  an  approximately  uniformly  distributed  0-1  matrix  with 
fixed  row  and  column  sums  is  by  defining  an  irreducible  Markov  chain  {A/,}  on  the 
interchange  graph  with  uniform  stationary  distribution.  If  T  is  a  function  defined 
on  0-1  matrices,  then 

1  ” 

-  y  r(A/.)  ET 

11  ' 
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with  probability  1.  If  Mo  is  a  fixed  matrix  in  this  set,  we  can  then  estimate  P[T  > 
T(Mo))  by  simulating  the  Markov  chain  and  computing  empirical  probabilities  after 
a  suitably  large  number  of  transitions.  In  particular,  let  D  be  the  meiximum  degree 
of  any  vertex  in  the  interchange  graph,  let  p  <  D~'^^  and  for  a  vertex  M,  define 
P(M,M')=p  if  {M,M^}  is  an  edge 
P(M,M)  =  l~Dp 
P(M,M')  =  0  otherwise 

It  may  be  shown  that  these  transitions  define  an  irreducible  Markov  chain  on  the 
interchange  graph  with  uniform  stationary  distribution.  Call  such  a  chain  a  modified 
random  walk. 

The  effectiveness  of  this  approach  depends  on  how  rapidly  the  distribution  of  M,- 
converges  to  the  uniform  distribution.  For  binary  matrices  with  fixed  row  and  col¬ 
umn  sums,  the  state  space  may  be  very  large  indeed.  In  general,  it  is  approximately 
a  hyper-exponential  function  of  the  size  of  the  matrix  M.  See  Good  and  Crook  [4] 
for  asymptotic  approximations,  as  well  as  exact  recurrence  formulas  for  certain  spe¬ 
cial  cases.  Ideally,  we  want  the  number  of  steps  needed  for  the  distribution  of  M,- 
to  converge  to  be  a  polynomial  in  the  “size”  of  M.  The  main  result  of  this  paper 
is  that  for  a  pai  ticular  modified  random  walk  on  the  interchange  graph,  the  eigen 
value  defining  this  “relaxation  time”  is  bounded  by  a  polynomial  in  the  number  of 
non-zero  entries  and  the  number  of  rows  and  columns  in  M,  but  is  otherwise  inde¬ 
pendent  of  the  number  of  0-1  matrices  satisfying  the  row-column  sum  restrictions. 

Binary  matrices  with  row  and  columns  sums  fixed  can  be  regarded  as  matrix  rep¬ 
resentations  of  bipartite  graphs  with  a  given  degree  sequence.  The  problem  of  sim¬ 
ulating  general  graphs  with  a  fixed  degree  sequence  can  be  approached  by  similar 
methods.  Here,  one  constructs  a  Markov  chain  on  a  space  of  graphs  by  selectively 
adding,  removing,  or  exchanging  individual  edges  in  the  graph.  Again,  under  cer¬ 
tain  restrictions  on  the  degree  sequence,  we  show  that  the  subdominant  eigenv^due 
of  this  chain  is  bounded  by  a  polynomial  in  the  total  number  of  edges. 

Uihforin  distributions  on  spaces  of  integer  matrices  subject  to  vTirious  row  and 
column  constraints  arise  in  a  number  of  areas.  A  very  specific  example  is  in  biogeo- 
graj)hiciil  ecology.  Here  we  have  a  collection  of  habitats  under  study  and  a  collection 
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of  species  of  interest.  For  each  habitat,  we  know  which  species  occur  in  that  habitat. 
If  we  write  the  data  as  an  incidence  matrix  of  species  in  habitats,  we  get  a  binary 
matrix,  say  Mq. 

Ecologists  are  particularly  interested  in  whether  or  not  two  species  compete.  One 
way  of  assessing  the  degree  of  competition  between  two  species  i  and  j  is  to  count 
the  number  tij  of  habitats  where  both  species  occur  and  then  compute  the  prob¬ 
ability  that  a  random  arrcingement  of  species  in  habitats  would  have  tij  or  fewer 
common  habitats  for  species  i  and  j.  We  restrict  the  random  arrgmgements  by  re¬ 
quiring  that  in  any  such  arrangement,  each  species  must  occupy  the  same  number 
of  habitats  and  each  habitat  must  support  the  same  number  of  species  as  in  our 
original  Mq.  Then  the  desired  probability  is  P[Tij{M)  >  where  M  is  uniformly 
distributed  over  binary  matrices  with  row  and  column  sums  equal  to  those  of  Mq. 
See  Simberloff[9],  Conner  and  Simberloff[3],  and  Simberloff  and  Zaman  [10]  for  more 
detailed  discussion  of  this  problem. 

Another  place  where  random  integer  matrices  arise  is  in  the  study  of  contingency  ta¬ 
bles.  Here,  conventional  statistical  analysis  might  begin  by  testing  for  independence 
of  rows  and  columns  by  means  of  the  statistic.  If  the  hypothesis  of  independence 
is  strongly  rejected,  there  is  a  need  for  some  alternative  probability  model  to  de¬ 
scribe  the  table. 

One  alternative  that  has  been  proposed  is  a  uniform  distribution  on  the  space  of 
contingency  tables  with  a  given  set  of  marginal  totals.  In  the  setting  of  the  present 
paper,  the  state  space  may  be  regarded  as  the  set  of  contingency  tables  with  the 
further  restriction  that  all  entries  be  either  0  or  1.  Diaconis  and  Efron  [5]  discuss  this 
model  for  contingency  tables,  along  with  a  family  of  other  models.  They  provide 
a  formula  for  computing  approximate  probabilities  for  this  distribution.  Markov 
cliain  simulations  using  the  interchange  formula  provide  an  alternative  means  for 
computing  probabilities  for  this  distribution. 

In  Section  2  of  this  paper,  we  review  some  standard  definitions  and  notation  for 
binary  matrices,  graph  theory,  and  Markov  chains.  In  Section  3,  wc  prove  an 
essential  graph  theory  lemma,  and  then  apply  it  to  bounding  the  cigcn\'alue  of  a 
Markov  chain  for  generating  an  almost  uniformly  distrib\iied  graph  with  a  fixed 
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degree  sequence.  The  main  results  of  this  paper  are  Propositions  6  and  7  in  section 
4,  where  we  compute  a  bound  for  the  rate  of  convergence  of  a  random  walk  on 
the  interchange  graph  for  a  space  of  binary  matrices.  Afterwards^  we  discuss  the 
convergence  rates  of  functions  of  the  chain. 


2.  Some  Definitions: 

For  integer  sequences  r  =  (ri,...  ,r„,)  and  s  =  (si,.. .  ,<Sn)i  let  be  2l(r,s)  be  the 
set  of  binary  matrices  with  row  sums  r  and  column  sums  s,  and  suppose  that 
2l(r,s)  0. 

Define  a  graph  J  =  (21,  S)^  with  21  =  ^(r,  s),  and  edge  set  C  defined  by  {M,  M'}  6  ^ 
if  there  exist  rows  ii  and  12  and  columns  ji  and  72  such  that 


Af[ii,t2;  ji,  j2] 


M[iui2]juj2] 


1  0’ 

1  f  *  *  1 

0  l" 

0  1 

M  [ii,i2;;i,;2l  = 

1  0 

0  1] 

%  t  •  *  1 

1  0 

— 

1  Oi 

L  4 

M  [ii,i2;;i,;2]  = 

0  1 

or 


(1) 


and  M[kJ]  =  M‘[kJ]^  for  k  ^  11,12  or  I  ^  jij2-  Call  the  operation  taking  M  to 
an  interchange,  and  J  the  interchange  graph  on  2J(r,  s).  We  will  often  write 
either  of  these  two  interchanges  in  the  form  ^2). 

One  important  interpretation  of  2l(r,s)  is  as  the  set  of  bipartite  graphs  with  biparti- 
tion  {!,...  ,m;  m  +  l, . . . ,  m-fn}  and  degree  sequence  ri, . . .  ,r^;  si , . . .  ,Sn.  We  shall 
use  either  the  graph  or  the  matrix  interpretation  of  21,  according  to  convenience. 

For  binary  matrices  A  =  [a,j]  and  B  =  [6,;],  let 


.4  0  J?  =  (a,j  4-  (mod  2)]  .4  V  B  =  (max{a,j, 6,^}  ]. 


For  future  reference,  we  will  review  some  standard  graph  theory  notions.  Let  G  = 
{V,S)  be  an  arbitrary  graph.  Say  that  7  =  fxo,xi,. . .  ,x*,rk4.i  =  xo}  is  a  circuit 
if  Xo, . . . ,  X*  arc  vertices  and  {x,,  x,>i }  is  an  edge  for  i  =  0, . . . ,  k.  If,  in  addition, 
Xo, . . .  ,Xjt  arc  distinct  we  say  that  7  is  eiemeiitary.  Say  that  G  is  Eulerian  if  there 
exists  a  circuit  0  that  traverses  every  edge  in  S  exactly  once;  we  call  the  circuit  ^ 
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an  Eulerian  circuit.  We  recall  the  elementary  result  that  a  graph  is  Eulerian  if  and 
only  if  it  is  connected  and  contains  either  no  or  two  vertices  of  odd  degree. 

The  fundamental  result  on  interchanges  is  that  the  interchange  graph  is  connected, 
known  as  Ryser’s  Theorem. 

Theorem  1.  Let  A  and  B  be  matrices  in  2l{r,s)  and  let  B— A  =  Ci+ — hC^,  where 
Cl, . . . ,  Cg  are  disjoint  elementary  circuits.  Let  the  number  of  non-zero  entries  of 

Ci  be  2ki,  i  =  1, . . . ,  g.  Then  there  exists  a  sequence  ofki^i - hk^—g  interchanges 

which  transforms  -4  into  L. 


Proof:  This  is  Theorem  3.2  in  Brualdi[2]. 

We  will  also  need  some  standard  ideas  about  Markov  chains.  Let  P  be  an  irreducible 
aperiodic  Markov  transition  matrix  on  some  finite  set.  Then  it  is  well-known  that 
P  has  a  unique  stationary  distribution  tt,  and  ||P”(x,o)  -  Tr(o)||  <  where  is 
the  subdominant  eigenvalue  of  P. 

Given  some  knowledge  of  the  geometry  of  our  state  space,  we  can  bound  /?i.  Let 
P  be  reversible,  so  that  7r(i)P(x,y)  =  7r(y)P(y,x)  for  all  states  x  and  y.  Let 
Q(x,y)  =  7r(x)P(x,y).  For  each  pair  {x,y}  in  the  state  space,  suppose  there  is  a 
path  7ry  connecting  x  to  y.  Let 

r/ =  maxQ(e)”'  ^  ;r(i)jr(y). 

Then,  the  following  theorem  holds. 


Theorem  2.  For  a  reversible,  irreducible  Afarkov  cliain  P,  the  second  largest  eigen 


value  satisfies 


<  1 


Proof:  This  is  Proposition  7  in  Diaconis  and  Stroock[Cj. 
Finally,  for  arbitrary  matrices  A,B,  define  the  matrix 


If  Ai , . . . ,  A„  are  real  numbers,  diag(Ai , . . . ,  A„)  will  denote  the  n  x  n  diagonal  matrix 
with  entries  Ai, . . . ,  A„. 

3.  Generating  Graphs  by  Edge  Perturbations: 

Let  d  =  (di, . . .  ,dn)  be  a  sequence  of  positive  integers.  Let  5(d)  =  5  be  the  set  of 
graphs  with  degree  sequence  d,  and  suppose  that  5^0-  Let  Q'  be  the  set  of  graphs 
with  degree  sequences  d'  =  (dj, . . .  ,d^^)  satisfying 
X.  d\  <  di,  i  =  1, . . .  ,n 
n.  -d\)  =  2 

Let  5  =  5  U  5^  In  [8],  Section  2,  Jerriim  and  Sinclair  define  the  following  set  of 
transitions  for  a  Markov  chain  on  S. 

i.  Select  an  edge  {z,i}  uniformly  at  random. 

:t.  U  G  £  Q  and  {i,i}  is  an  edge  in  G,  let  =  G  — 

Hi.  If  G  €  5^  {hj]  is  not  an  edge  in  G,  and  the  degree  of  i  is  less  than  d,,  let 
H  =  G 4-  {t, j }.  If  the  degree  of  j  exceeds  dj  select  an  edge  {j, k]  uniformly  at 
random  and  delete  it. 

XV.  In  all  other  cases,  do  nothing. 

Let  d„,„jp  =  max(di , . . . ,  d„),  and  let  AT  =  “  dj  be  the  number  of  edges  in  G.  Say 
that  d  s  p’Stable  if  N  >  -  d,„ax.  In  (8),  Jerrum  and  Sinclair  remark  without 

proof  tiiat  a  Mzirkov  transition  matrix  defined  by  i.  -  xv.  will  have  a  subdominant 
eigenvalue  that  is  bounded  by  a  polynoniiiil  in  (dj , . . . ,  d„)  if  d  is  p-stablc.  We  now 
show  this  explicitly. 

Proposition  3.  The  subdoinimuit  eigcjxvaJue  of  the  transition  matrix  induced  by 
i.-iv.  satisfies  the  bound 

d,  <  1  -  i 

Proof:  Our  approach  follows  that  of  Jerrum  and  Sinclair  [7],  [8].  As  in  these 
papers,  we  will  define  a  system  of  canonicitl  paths  in  5,  use  these  to  bound  ij  and 
then  estimate  the  subdominant  eigeiu'alue  of  the  transition  matrix.  As  a  first  step 
we  prove  a  graph  theory  Icinina. 

Lcnitna  4.  Let  G  -  (1  be  a  connected  graph  Suppose  S  =  Si  U  S2,  Si  O 


G 


S2  =  0,  and  that  |ei(t;)|  =  lc2(v)|  for  aii  v  €  V,  where  ei(v)  and  €2(v)  are  the 
edges  in  Si  and  82^  respectively^  incident  to  v.  Then  G  has  an  Eulerian  circuit 
{vo,Vi,...  ,V2n-^i,vo}  such  that  {v2i,V2i+i}  e  Si  and  {v2i+i,  1^21+2)  €  S2  fori  = 
0,...,n  -  1. 

Proof:  It  is  easy  to  see  that  |£^|  is  even,  say  \S\  =  2n.  The  proof  follows  by  an 
induction  on  n.  The  smallest  n  for  which  the  conditions  of  the  lemma  can  be 
satisfied  is  n  =  2.  For  n  =  2,  we  have  V  =  {1,2, 3, 4},  Si  =  {{1, 2},  (S,  4}}  and 
S2  =  {{2,3},  {4, 1}},  within  isomorphism.  7  =  {1,2, 3,4, 1}  is  then  the  desired 
circuit. 

Now  suppose  the  lemma  holds  for  n  =  Jb,  and  let  n  =  ib  -h  1,  so  |£|  =  2^  -f  2.  Choose 
a  vertex  t’o  and  edges  {uo,t;i}  and  {ui,t;2}  in  Si  and  82^  respectively.  These  exist, 
by  hypothesis.  Continuing  in  the  same  manner,  select  vertices  U3, . . .  ,i;2m  so  that 
{t^2i,  i’2i+i }  €  {t^2i>ii  1^21+2)  €  £^2»  for  i  =  1, . . . ,  m  -  1  and  all  edges  are  distinct, 
until  you  reach  the  first  vertex  V2m  where  such  a  choice  cannot  be  made.  I  claim 
for  if  {u2m,u^}  appears  in  the  circuit  for  all  w  adjacent  to  U2m,  then  an 
odd  number  of  edges  adjacdit  to  V2m  must  have  been  crossed  prior  to  stage  2m. 
The  only  vertex  in  the  grapa  satisfying  this  requirement  is  t^o- 

Let  7  =  {eo,ri, . . .  ,i’2m-i,t^2m}*  If  m  =  fc  +  1,  then  7  is  the  desired  circuit,  zmd 
we  are  finished.  So,  suppose  m  <  ib  f  1.  Let  T  =  {{r,-i,t^i},*  =  i,...,2m},  let 
S^  —  S\T  and  let  G'  be  the  subgraph  induced  by  S\  Write  G'  =  Gj  U  •  •  •  U  G*, 
as  the  union  of  connected  components.  Then  it  is  easy  to  see  that  each  component 
G,  satisfies  the  hypotheses  of  this  lemma,  with  fewer  than  2ib  edges.  Thus,  suitable 
cir':uits  can  be  coi..  tructed  on  each  G,  and  adjoined  to  7.  The  lemma  follows.  ■ 

Corollary  5.  Let  G  =  (V,S)  he  a  connected  graph  as  in  Lemma  4.  Suppose  there 
exist  vertices  t  V'  such  that  |ci(r+)|  =  |f3(r^)|-f  1  aiid  )|  =  |c2(i’-.)|- 

1;  for  all  other  vertices,  suppose  =  €2(1*).  Then,  G  has  an  Eulerian  path 

{t‘o,  t’j , . . .  ,  f„ }  such  that  I’o  =  r„  =  r_  ,  ojid  {t’21,  i’2i+i }  t  Li,  {e2i+i ,  V2,+  2}  € 
S2  for  i  =  0, ...  ,Ti  -  1. 

We  now  proceed  lo  construct  our  set  F  of  canonical  paths.  Let  G,  //  6  C  and  let 
D  =  G  H  be  the  symmetric  difference  of  G  and  H .  Write  D  =  Di  U  •  •  V  Dt 
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where  Di , . . . ,  Djt  are  the  connected  components  of  D.  Any  vertex  v  in  D  has  the 
property  that  the  numbers  of  edges  in  G\H  and  H  \  G  incident  to  v  are  the  same. 
Thus  Lemma  4  applies  to  each  component  D,*. 

Order  the  Eulerian  subgraphs  of  G  U  in  some  fashion,  and  for  each  such  Eulerian 
subgraph,  suppose  a  starting  vertex  uq  has  been  identihed.  Suppose  that  Di , . . . ,  Dik 
is  an  increasing  sequence  in  the  order. 

The  canonical  path  from  G  to  H  will  be  defined  by  unwinding  , . . . ,  Dik  in  order. 
The  unwinding  will  be  carried  out  as  follows: 

>.  For  each  component  £),,  fix  a  starting  vertex  v.  Let  X>i  be  an  ordering  of  the 
Eulerian  circuits  on  D,,  with  v  as  starting  (and  ending)  vertex.  Note  that 
every  Eulerian  circuit  induces  an  ordering  on  the  edges  of  E,. 
it.  Let  ® ,  be  the  set  of  Eulerian  circuits  on  Di  such  that  edges  in  G\  If  have  odd 
parity  iuid  edges  in  H  \  G  have  even  parity.  (That  is,  if  7  =  {t'o,*  •  • 
then  {t;2,, U2t-n}  ^  G  \  H  and  {v2i-f i,U2,>2}  G  /f  \  G  for  i  from  0  to  m  —  1.) 
By  Lemma  4,  is  not  empty.  Let  /?  be  the  first  element  in  ©,,  with  respect 
to  our  ordering  of  Pj. 

hi.  Suppose  p  =  {vo,. . .  ,V2m}.  The  unwinding  begins  by  removing  {uo,iq}.  The 
next  m  ~  1  steps  consist  of  adding  {v2i~i,*^2i}  ^nd  lemoving  {v2nV2i-fi)  in 
sequence  for  i  =  1, . . . ,  rn  -  1.  The  finiil  step  is  to  add  the  edge  {v2m~i)i^2m}- 
Let  be  a  transition  in  the  cemonical  path  from  G  to  H,  and  let  v  be  a 

vertex  in  M  0  M' .  By  the  definition  of  canonical  paths,  cg(v)  fl  c/f(v)  C  CA/nAf'(v) 
C  C(:;(e)  U  c Let  5  =  A/  U  A/'  0  (G  0  H).  By  elementary  set  theory, 

k.v(r’)l  =  kA/uAf'(v)0(c(7(i')0c//(e))| 

=  kA/u.\/'(t’)  U  (C(;(r)  0  c//(v))|  ~  |cAfuAr'(t^)  Fl  (tc;(r)  c//(v))| 

=  |rc;(i’)U  r//(e)l  -  (k.\foAf'(v)|  -  |cf;(r)  0  c//(r)|) 

=  k(;(i')l  +  \ct{(v)\  ~  k.«u.\f'(v)| 

There  are  three  possible  cases  to  consitler.  Let  i  be  the  unicpic  vertex  of  degree 
2  in  M  0  Then  ir .\fnA/'(i’)l  -  d,  -f  1,  and  as  G, //  G  Qy  k5(i))|  =  ti,  -  1. 
Alternatively,  let  j  be  the  starting  vertex  in  the  cycle  7  we  are  unwinding.  Then 
kAf(;  )l  -  kA/'{;  )l  -  -  1,  so  k-v{;))|  -  d^  -f  1.  For  all  other  vertices  k,  ko'(i’)|  = 

k//(^‘)l  -  kA/nAf'(i’)|  “  di,  so  k.s*(*'’))|  =  di. 
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Define  a  function  (7r{G^H)  encoding  G  and  H  by  setting 

TT\  _  j  {G  ^  H)  ®  (M  U  M')  if  r  begins  or  ends  a  circuit 

^  “  I  (G  0  if)  0  (M  U  M')  -  eG,<  otherwise 

Here,  eG,t  is  the  first  edge  in  13. 

From  the  discussion  in  the  last  paragraph,  it  follows  that  ar(G,if)  €  S.  Further¬ 
more,  a  is  one-to-one,  since  we  can  reconstruct  G  and  H  from  <Tr(Cj,  H)  amd  AiUAf 
To  see  this,  let  K  —  ar{G,H)^(M UM').  If  ar{G^H)  €  G,  then  K  =  Di  U  -  •  -  UD*. 
If  (7r(G,  H)eg\  then  iv  =  Dj  U  •  •  •  U  A-i  U  Di  U  D.+i  U  •  •  •  U  Dit,  where  Di  has 
precisely  two  vertices  of  odd  degree.  We  recover  Di  by  adding  an  edge  joining  the 
two  odd-degree  vertices.  D,  is  the  unique  component  containing  M  fl  so  D,  is 
the  current  subgraph  in  the  unwinding. 

Let  C,  be  the  set  of  Eulcrian  paths  on  D,  starting  at  the  initial  vertex  v,  and  having 
the  property  that  edges  in  <Tr(G,  H)  have  odd  parity  prior  to  M  0  M'  and  even 
parity  afterwards.  /?  €  C,  so  C,  is  not  empty.  Let  ^  be  the  first  path  in  G,.  Then 
[3  induces  a  path  on  D,  such  that  edges  in  G  and  H  have  odd  and  even  parity 
respectively.  Thus,  0  G  and  H  can  now  be  identified  by  unwinding  0  forwards 
or  backwards  from  M.  as  necessary. 

To  extend  canonical  paths  to  the  whole  of  5,  we  will  define  a  function  G  G  that 
will  associate  a  “nearest"'  member  of  Q  to  each  clement  in  5.  If  G  €  let  (5  =  G. 
Now,  let  G  ^  G'.  Following  Jerrum  and  Sinclair [S],  define  a  graph  (5  €  5  that  is 
“close”  to  G.  There  are  two  cases  to  consider. 

i.  Suppose  there  are  vertices  i  and  j  such  that  =  cf,  —  1  and  =  dj  —  1.  If 
{*»)}  i  (5  =  G  -f  {Lj}.  If  ^  S,  then  find  an  edge  {kj}  such  that 
{z,  It}  iukI  {jj]  are  not  edges  and  let  (5  =  G  -f  {i,  k)  -f-  {;,  /}  -  {^‘,0-  shown 
in  Jernnn  aiul  Sinclair[S],  section  3,  such  an  edge  zdways  exists,  if  the  sequence 
d  is  p  stable. 

ii.  Alternatively,  suppose  there  is  a  vertex  i  such  that  =  d,  -  2.  Again,  find  zui 
edge  [kj]  such  that  {i,  I‘)  and  {ij]  arc  not  edges  and  let  G  =  G  -f  {t,  d* 
{ul]~{kj}. 

In  either  case,  since  G  and  G  differ  by  at  most  two  edges,  at  most  A  *  graphs  G  can 
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be  associated  with  any  G. 

For  Gy  H  €  Sy  define  the  canonical  path  from  G  to  if  by  connecting  G  to  G,  if  to 
Hy  and  defining  the  canonical  path  from  G  to  J?  as  in  the  preceding  paragraphs. 
Then,  the  number  of  canonical  paths  crossing  any  given  edge  {M,  M'}  is  at  most 
N^\G’\.  As  before,  let  t]  =  maXcQ(e)"^ 

V  =  Q(e)->  max  E  i4  ^  2Af|5|  x  N*\S\  x =  2N^ 

7*»3«  '  *  '  ' 


Apply  Theorem  2  to  get 

This  completes  the  proof  of  the  proposition.  ■ 


Remark:  In  [8],  Jerrum  and  Sinclair  describe  an  alternative  algorithm  for  almost 
generating  uniform  random  graphs  in  5(d)  by  tranlating  the  problem  into  one 
of  simulating  random  perfect  matchings.  Additional  structure  provided  by  the 
matching  problem  gives  the  bound 


Pi<l- 


_ i _ fMVV 

16n''(n/2  +  2Af)(n-l)  Vie'iy  J 


for  an  arbitrary  family  5(cl).  The  p-stability  condition  gives  a  sufficient  condition 
for  151/151  to  be  polynomially  bounded.  If  d  is  p-stable,  then 

/  1 

/?i  <1- 


^16iV4n4(n/2-f  2A^)(n-l)^ 


4.  Simulating  Binary  Matrices: 

Let  r  =  (ri , . . . ,  7-„, )  and  s  =  , . . . ,  )  be  integer  sequences  such  that  2l(r,  s)  is 

not  empty.  For  convenience,  write  21  =  2l(r,s). 

VVe  will  need  notation  for  some  other  spaces  of  binary  matrices.  Let  21^  =  U2l(r1s^), 
where  the  union  is  taken  over  sequences  r'  =  (r'j, . . . ,  s'  =  (sj , . . . ,  s'„)  such  that 

=  r,  -  1  rj.  =  rjk, 

I  y  t 

Sj  =5^-1  S,  =  Si, 
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for  some  indices  i  and  j.  Similarly,  let  21^^  =  U2l(r”,s"),  where  r"  =  (r”,. . . 
s'' =  (3'/,..., s''),  and 


for  pairs  of  indices  t'lj  ^2  stnd  ii,i2-  Let  5  =  21 U  21^  U  21". 

Let  M  G  21.  Define  transitions  for  a  Markov  chedn  on  21  as  follows: 

i.  Select  two  entries  m^v  and  rurs  uniformly  at  random  from  the  entries  in  M 
equal  to  1. 

u.  If  771,^3  =  TTirv  =  0,  interchange  rriuv  a**  i  rrirs  with  rrius  and  nirv 
in.  Otherwise,  do  nothing. 

Let  P  be  a  transition  matrix  corresponding  to  this  transition  scheme.  Let  7r(M)  = 
M  G  21.  Then,  it  is  not  hard  to  show  that  tt  satisfies  the  detailed  balance 
equations  for  P.  If  {Xn}  is  a  Markov  chain  on  21  with  these  transitions  then  {Xn} 
is  reversible,  with  stationary  measure  tt. 

We  wish  to  estimate  the  rate  of  convergence  of  the  distribution  of  Xn  to  tt.  As  in 
Section  3,  we  will  construct  a  system  of  canonical  paths  in  21  .  The  construction 
will  be  done  so  that  for  any  matrices  J  and  P  in  21  and  any  edge  {M,  M'}  lying  on 
the  canonical  path  from  J  to  P,  we  can  map  (J,  P)  uniquely  to  a  matrix  in  S. 

Proposition  6.  There  exists  a  set  V  of  paths  joining  every  pair  of  points  in 
21  such  that  the  number  of  paths  containing  a  hxed  edge  is  bounded  above  by 

2(("V')  +  l)|5i. 

Proof:  Let  J,  P  G  21  and  suppose  J  ^  F  =  C,  where  C  is  an  elementary  circuit. 
There  is  a  standard  algorithm  for  transforming  J  to  F  by  intercheinges.  (See, 
for  example,  Bnialdi[2],  Section  3.)  Let  C  have  edges  {xj , yi },  {yi , i2}>  y2}» 

. . . ,  {xt,  yfc },  {t/A:,  xi }.  Without  loss  of  generality,  suppose  the  matrices  representing 
J  and  P  have  Jjj  =  J22  =  •••  =  Jkk  =  I»  J12  ~  *^23  =  —  Jk-\,k  =  ^ki  — 

Fn  =  P22  =  ...  =  Fkk  =  0,  Pi2  =  P23  =  • . .  =  Fk-i,k  =  Pi-i  =  1,  and  =  P,j 
elsewhere.  Let  p  be  the  smallest  integer  such  that  Jpi  =  0;  Jk\  ==  0,  so  p  <  k. 
The  first  p  interchanges  will  be  (p,p  -  1;  l,p),(p  -  l,p  -  2;  l,p  -  1), ...  (2, 1;  1,2). 
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Figure  1.  The  matrices  (a)  J,  (b)  F,  (c)  J  ©  F,  and  (d)  J\  where  m  =  n  =  7, 
p  =  4 

These  p  interchanges  produce  a  matrix  J'  where  =  ‘/p+2,p+2  =  . . .  = 

Jp+i^p^2  =  ...  =  ‘^I'j  “  ^ij  elsewhere.  (See  Figure  1  for 

typical  examples  of  J,  F,  J  ©  F,  and  J^)  Repeating  this  process  will  eventually 
transform  J  to  F.  Note  that  this  algorithm  proceeds  by  unwinding  a  series  of 
elementary  circuits  which  are  contained  in  C,  with  the  possible  exception  of  one 
edge.  Call  these  circuits  subsidiary. 

It  will  be  convenient  to  represent  canonical  paths  somewhat  differently.  Associate 
each  matrix  .4  6  21  with  the  matrix  A  =  diag(M,  1,1).  We  can  map  the  canonical 
path  from  J  to  F  to  a  somewhat  different  path  from  J  to  F  As  before,  let  p  be 
the  smallest  integer  such  that  Jpi  =  0.  Let  the  first  interchange  in  the  canonical 
path  from  J  to  F  be  (l,m  -1-  1;  l,n  -1-  1).  If  p  =  2,  let  the  second  interchange  be 
(l,2;n  +  1,2);  otherwise,  let  the  second  interchange  be  (p,  rn-f  2;p,  n4-2).  The  next 
p-*  1  interchfuiges  will  be  (m +  2,p- l;p,p~  l),(m>f2,p~2;p-’  l,p~2),.. .  ,(.714- 
2,1;2,1),(1  ,  m  -f  2;  l,n  4-  1);  these  correspond  precisely  to  the  interchanges  (p,p 
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1;  l,p), (p  —  l,p  —  2;  l,p*-  1), . . . , (2, 1;  1,2).  Finally,  make  the  interchange  (l,p; n -f- 
l,n  -f  2).  The  unwinding  will  continue  using  the  (p,  n  -f  1)  entry. 


Let  be  an  edge  in  the  canonical  path  from  J  to  F.  Without  loss  of 

generality,  write  the  corresponding  interchange  as  (/,/  —  translate 

this  to  an  interchange  on  the  modified  canonical  path  in  a  standard  fashion.  Suppose 
{M,  M'}  lies  on  a  subsidiary  cycle,  with  initial  entry  q  and  final  entry  p.  Transform 
M  and  M'  by  the  interchanges  l,n-f  1),  (/,  m-i-2; /, n  +  2),  (/,p-hl;  l,n4-l) 

and  1;  l,n  + 1),  (/ —  1,  m -f  2; /  —  l,n -f- 2),  (/ —  l,p-i- 1;  1,  n -f  1),  respectively. 

This  will  transform  M  and  M'  into  matrices  M*  and  M*'  such  that  M*  is  joined  to 
M*'  by  the  interchange  (/  —  1,  m  +  2;  /,  /  —  1).  Alternatively,  suppose  {M,  M'}  does 
not  lie  on  a  subsidiary  cycle.  Then  it  is  easy  to  see  that  (/,  /  —  1;  1,  /)  corresponds 
to  the  interchange  (/,/  —  l;n  +  1,/).  To  summarize,  each  edge  {M, M'}  in  the 
interchange  graph  can  be  transformed  to  a  modified  edge  by  either  one 

or  three  interchanges. 


The  first  p  -f-  1  interchanges  transform  J  into  a  matrix  where  Jp^2,p+2  =  •  •  •  = 

^kk  ^  •^p+i,p+2  ~  •^p-f2,pT3  =  . . .  =  ^k~i,k  ^ki  —  ~  other 

l<z<m,  l<.7<n.  Again,  we  can  repeat  this  process  to  transform  J  into  F 
eventually. 


Now,  suppose  J  0  F  has  connected  components  Gi,. . .  ,Gr,  and  let  G  be  one  of 
these  components.  As  in  Section  3,  specify  a  starting  vertex  for  each  G,-  and  order 
the  Eiilerian  circuits  starting  from  that  vertex.  By  Lemma  4,  each  G,  has  at  least 
one  Eulerian  circuit  in  which  entries  from  J  and  F  alternate.  Let  7  be  the  first 
such  circuit.  Write  7  as  . .  ,Xj,,y^,ii}.  Let  p  =  inin{j  :  Xj  =  i,, t  < 

j}  A  inin{j  :  yj  =  y,, i  <  j},  and  let  j{p)  be  the  unique  j  <  p  such  that  Xj  =  Xp 
or  yj  -  y^.  Then  ,  (yp_,,x^}  or  {yi(/,),Xj(p)+i  },•••,  as 

the  ciisc  may  be,  forms  an  elementary  circuit,  which  we  denote  Ci.  Repeating  this 
argument  successively  on  G  \  Ci,  G  \  (Ci  U  C2)  and  so  forth,  we  come  to  an  ordered 
decomposition  of  G  into  elementiuy  circuits  Ci, . . . ,  Ck- 

Since  7  is  an  alternating  Eulerian  circuit,  each  C,  is  an  alternating  elementary  cir> 
cult.  We  can  piiss  from  J  to  F  by  successively  unwinding  Gj, . . . ,  Gr-  In  turn,  each 
component  G j  can  be  unwound  by  unwinding  the  elementary  circuits  , . . . ,  Cj  k, 
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in  order.  So  the  unwinding  algorithm  for  elementary  circuitf.  extends  to  an  unwind¬ 
ing  algorithm  for  Eulerian  circuits,  thence  for  the  difference  between  elements  in 
Ql. 

To  count  canonical  paths  containing  we  will  first  encode  modified  paths 

containing  as  a  pair  (p,cr)  G  {+1,-1}  x  5.  Let  J  =  diag(J,  1,1),  F  — 

diag(F,  1,1),  and  let  5(J,  F)  be  the  principd  m  x  n  submatrix  of  J©F©(M*' VM*). 
The  degree  sequence  of  5  will  depend  on  the  position  of  {M*^ ,  M* }  in  the  unwinding. 
There  are  three  possible  cases: 

i.  ends  zin  elementary  circuit. 

ii.  {M*\M*]  lies  inside  an  elementary  circuit,  but  is  not  on  a  subsidiary  circuit. 
Let  X,  be  the  starting  vertex  of  the  elementary  circuit,  and  let  yj  be  the 
current  vertex  being  changed.  Calculations  similar  to  those  in  section  3  show 
that  deg5(x,)  =  r,-  —  1,  deg5(yj)  =  sj  +  1.  Furthermore,  it  is  not  hard  to  see 
that  the  first  edge  (2,2^  of  the  current  elementary  circuit  is  in  S. 

Hi.  lies  on  a  subsidiary  circuit,  with  starting  edge  {p,p'}  and  ending 

As  before,  we  will  have  deg5(x,)  =  r,-  —  1,  degsiVj)  —  +  1- 

Also  deg5(xp)  =  Tp  +  1,  and  deg^{xq)  =  r,  +  1,  and  the  edges  {p,p'},  and 
{q,q^}  are  in  5. 

With  this  in  mind,  define 


cr{j\F)‘-= 


<  S(J,F) 

l5(J,F) 


€i,i' 

Ci,.' 


m  case  1. 
in  case  ii. 
in  case  Hi 


a!iO 


I  ri  lu  cases  t.  and  ii. 

p{J ,F)  =  <  +1  in  case  in’.,  if  (p  -  q){p'  -  q‘)  >0 
{  -1  in  case  iii,,  if  (p  -  q){p‘  -  q‘)  <  0. 


Here,  c,,,s  Cpy,  and  Cq  q>  are  the  edges  joining  {i,i^},  {p,p^},  and  {q^q']^  respec¬ 
tively. 

Exmnining  the  definition  will  show  that  a{J,  F)  G  21  in  case  i.,  a{J,  F)  G  21'  in  case 
ii.,  and  a(J,F)  G  2l"  in  case  iii. 
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We  claim  thc  t  we  can  determine  J  and  F  from  (/>,cr)  and  {M*M*^}.  To  see  this, 
first  observe  that  we  can  reconstruct  S{J,F)  from  (/),  <t).  In  cases  i.  and  n.,  this  is 
identical  to  the  argument  in  Section  3.  In  case  zu‘.,  we  have  to  identify  the  edges 
Cpy,  and  Cq^qi  and  restore  them  to  cr  to  obtain  5.. 

Since  there  are  precisely  two  rows  in  a  with  diminished  row  sums,  one  of  which 
corresponds  to  M*  0  M*\  and  since  M*[n  +  l,zM  =  1,  reconstructing  is  auto¬ 
matic.  To  reconstruct  Cpy  and  eg,g/,  note  that  M*[p,  n  +  1]  =  M*[g,  n  -1-2]  =  1, 
Sp,  =  Spi  —  1,  and  ^  Sqi  —  1.  This  leaves  two  possible  pairs  of  edges,  which  are 
distinguished  by  the  value  of  p. 

Knowing  S{J^F)^  we  can  reconstruct  J  and  F  by  an  argument  similar  to  that  in 
Section  3.  Let  J  0  F  have  connected  components  Gi, . . .  ,Gr.  There  is  a  unique 
component  containing  the  edges  {u,  u}  and  {u,  u'}.  For  j  ^  k,  let  jj  be  the  first 
Eulerian  path  on  Gj  such  that  edges  in  5  and  M*  V  M*'  alternate.  For  1  <  j  <  k, 
odd  parity  edges  are  in  J  and  even  parity  edges  are  in  F;  for  ac  <  j  <  r,  the  parity  of 
edges  is  reversed.  Let  be  the  first  Eulerian  path  on  Gk  such  that  its  elementary 
circuit  decomposition  Fi,...,Fp  gives  even  parity  to  edges  in  S  preceding  {u,u} 
and  odd  parity  to  edges  in  S  following  {v,u'}.  Odd  parity  edges  preceding  {u,t;} 
and  even  parity  edges  following  {u,  u'}  are  in  J;  Even  parity  edges  preceding  {u,  v} 
and  odd  parity  edges  following  {u,  u'}  are  in  F.  Thus,  J  and  F  can  be  reconstructed 
from  M*  V  M*'  by  adding  and  deleting  appropriate  edges. 

Since  (J,  F)  cr(J,F)  is  invertible,  it  follows  that  the  mapping  is  one-to-one. 
Each  encoding  cr(J,F)  is  in  S.  Thus,  the  number  of  canonical  paths  containing 
{M*\M*]  is  at  most  |5|. 

It  remains  to  compute  the  number  of  modified  transitions  {M*\M*)  corresponding 
to  each  possible  transition  {M',  A/}.  Depending  on  its  position  in  the  the  elemen¬ 
tary  circuit  F,  each  transition  {A/',  A/}  can  be  transformed  to  a  modified  transition 
{A/*',  A/* }  by  either  one  or  three  interchanges.  If  one  interchange  is  required,  then 
the  interchange  affects  one  of  the  two  vertices  in  A/\A/';  this  gives  two  possible  pairs 
{A/*', A/*}.  If  three  interchanges  are  required,  then  {M\M]  lies  on  a  subsidiary 
cycle.  Thus,  we  need  to  choose  the  beginning  and  the  end  of  the  subsidiary  cycle, 
by  choosing  any  two  ones  in  a  column  modified  by  {M\  A/}.  This  can  be  done  in  at 
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most  2  (^7^ )  possible  ways.  Thus  there  are  at  most  2  )  +  2  pairs  {M*',  M*} 

corresponding  to  each  pair  {M\M}.  The  total  number  of  canonical  paths  crossing 
{M\M]  is  then  at  most 


As  before,  the  estimate  of  the  number  of  canonical  paths  gives  an  estimate  of  the 
subdominant  eigenvalue  of  the  transition  matrix  of  {A’n}. 


Proposition  7. 


2/^  7  |21L 


-2 


Proof:  Again,  we  have  t]  =  maXe  Q(e)  ^  Using  our  bound  on 

the  number  of  canonical  paths. 

The  proposition  now  follows  from  Theorem  2.  ■ 


Proposition  8.  Let  21  =  2l(r,s).  Then  a  sufficient  condition  for  |5|/|2l|  to  be 
poIynoinidJy  bounded  is  that  N  >  2(r^ar  -  l)(5mox  -  1)  +  1. 

Proof:  Let  {r\ , . . . ,  r'„)  and  (sj , . . . ,  s'„)  be  the  row  and  column  sums  for  a  matrix 
hi  Cl'.  Then  -  1)  <  (rmar  -  iK^mar  -  l)-  SuppOSC  N  >  2{rm&x  - 

l){Srncr  -  1)  +  L  Then  N-l>  2{rmax  -  l)(^max  -  1)  >  2(r:„„  -  ^  1), 

SO  2l((r;,...,r;„);(s;,...,s;))  is  itself  a  p-stable  set.  As  the  number  of  p<  •  "ible 
.sequences  (rj , . . , ,  and  (sj , . . . ,  5',)  is  a  polynomial  in  m  and  n,  it  follows  that 
tlie  ]5|/|^3|  is  polynomially  bounded.  ■ 
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Our  estimate  of  /?i  gives  a  bound  on  the  average  error  committed  by  using  sample 
averages  of  our  Markov  chain  to  approximate  integrals.  We  state  our  bound  for  a 
general  chain. 

Let  {A'n}  be  a  Markov  chain  on  a  finite  state  space,  with  stationary  distribution  tt. 
Let  /  be  a  function  and  let  /i/  =  E/(2^)^r>  ~  For  m,n  >  0, 

let  Am^n  =  denotes  expectation  for  the  chain  with  Xq  =  x, 

then  the  following  estimate  is  routine: 


Proposition  9. 


"'a} 


(13) 


where  tt*  =  minj  7r(x). 


Proof:  For  m  ^  Poisson  (iV),  this  estimate  is  given  in  Aidous[l]  Proposition  4.2, 
with  in  place  of  and  some  differences  in  notation.  For  fixed  m,  it  suffices 
to  follow  the  argument  given  in  Aldous  to  his  expression  (4.9),  and  then  to  observe 
that  for  any  fixed  m  >  0  and  any  A,  A'”  =  (1  -f  (A  —  1))'"  <  exp(m(A  —  1)).  ■ 

Let  {A'n}  be  the  Markov  chain  on  21,  with  transition  matrix  P,  and  suppose  we 
wish  to  estimate  /i/  for  some  function  /;  here,  7r(x)  =  |2l|~L  To  estimate  /i/  with 
MSE  smaller  than  c,  we  can  take  m  >  log  M/ log /?i ,  giving  1  4*  |2l]/3{”  <  2.  Then 
take  n  >  8<7^/<(l  -  ^i)  to  give  Er{Am,n  -  ft <  e. 
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